<html>
<body style="width:800px;">
<h2>
    Profiling adventures and cython - setting the stage
</h2>

<p>
    <b>Summary</b>  This post is the first in a series dedicated to
    examining the use of profiling and, eventually, using cython, as a
    means to improve greatly the speed of an application.  The intended
    audience is for programmers who have never done any profiling and/or
    never used cython before. Note that we will not make use of cython
    until the third post in this series.
</p>


<img src="colored1.png" />

<h4>
    Preamble
</h4>

<p>
    Python is a great multi-purpose language which is really fun to use.
    However, it is sometimes too slow for some applications.  Since I only
    program for fun,
    I had never really faced a situation where I found
    Python's speed to be truly a limiting factor - at least, not
    until a few weeks ago when I did some exploration
    of a <a href="http://aroberge.blogspot.com/2009/12/17x17-challenge.html">
    four-colouring grid problem I talked about</a>.
    I started exploring ways to speed things up using only Python and trying
    to come up with different algorithms, but every one I tried was
    just too slow.  So, I decided it was time
    to take the plunge and do something different.  After considering
    various alternatives, like using
    <a href="http://code.google.com/p/shedskin/">shedskin</a> or
    attempting to write a C extension (which I found too daunting
    since I don't know C), I decided to try to use
    <a href="http://cython.org/">cython</a>.
</p>

<p>
    <a href="http://cython.org/">cython</a>, for those that don't know it,
    is a Python look-alike language that claims
    <i>
        to make writing C extensions for the Python language
        as easy as Python itself.
    </i>
    After looking at a few examples on the web, I concluded that such a rather bold
    statement might very well be true and that it was
    worthwhile trying it out on a more complex example.  Furthermore, I thought
    it might be of interest to record
    what I've done in a series of blog posts, as a more detailed
    example than what I had found so far on the web.
    As I was wondering if an esoteric problem like the
    four-colouring grid challenge mentioned previously was a good candidate
    to use as an example, by sheer serendipity, I came accross
    a <a href="http://www.reddit.com/r/Python/comments/ajeey/mandelbrot_set_in_python_tkinter/">link on reddit</a>
    by a new programmer about his
    <a href="http://prezjordan.tumblr.com/post/277984651/mandelbrot-set-in-python">simple Mandelbrot viewer</a>.
</p>

<p>
    Who does not like fractals? ... Furthermore,
    I have never written a fractal viewer.   This seemed like a good
    time to write one.  So, my goal at the end of this series of posts, is
    to have a "nice" (for some definition of "nice") fractal viewer that is
    fast enough for explorations of the
    <a href="http://en.wikipedia.org/wiki/Mandelbrot_set">Mandelbrot set</a>.
    In addition,
    in order to make it easy for anyone having Python on their system to
    follow along and try their own variation, I decided to stick by the
    following constraints:
</p>

<b>
<ul>
    <li>
        With the exception of cython, I will only use modules found in the
        standard library.   This means using Tkinter for the GUI.
    </li>
    <li>
        The program should work using Python 2.5+ (including Python 3).
    </li>
</ul>
</b>

<p>
    So, without further ado, and based on the example found on the
    reddit link I mentioned, here's a very basic fractal viewer that can
    be used as a starting point.
</p>

<pre>
''' mandel1.py

Mandelbrot set drawn in black and white.'''

import time

import sys
if sys.version_info < (3,):
    import Tkinter as tk
else:
    import tkinter as tk

def mandel(c):
    '''determines if a point is in the Mandelbrot set based on deciding if,
       after 20 iterations, the absolute value of the resulting number is
       greater or equal to 2.'''
    z = 0
    for iter in range(0, 20):
        z = z**2 + c
        if abs(z) >= 2:
            return False
    return abs(z) < 2

class Viewer(object):
    def __init__(self, parent, width=500, height=500):
        self.canvas = tk.Canvas(parent, width=width, height=height)
        self.width = width
        self.height = height

        # the following "shift" variables are used to center the drawing
        self.shift_x = 0.5*self.width
        self.shift_y = 0.5*self.height
        self.scale = 0.01

        self.canvas.pack()
        self.draw_fractal()

    def draw_pixel(self, x, y):
        '''Simulates drawing a given pixel in black by drawing a black line
           of length equal to one pixel.'''
        self.canvas.create_line(x, y, x+1, y, fill="black")

    def draw_fractal(self):
        '''draws a fractal picture'''
        begin = time.time()
        print("Inside draw_fractal.")
        for x in range(0, self.width):
            real = (x - self.shift_x)*self.scale
            for y in range(0, self.height):
                imag = (y - self.shift_y)*self.scale
                c = complex(real, imag)
                if mandel(c):
                    self.draw_pixel(x, self.height - y)
        print("Time taken for calculating and drawing = %s" %
                                                (time.time() - begin))

if __name__ == "__main__":
    print("Starting...")
    root = tk.Tk()
    app = Viewer(root)
    root.mainloop()
</pre>

<p>
    At this point, perhaps a few comments about the program might be useful
</p>

<ol>
    <li>
        I have tried to write the code in the most straightforward and
        pythonic way, with no thought given to making calculations fast.
        It should be remembered that this is just a starting point: first
        we make it work, then, if needed, we make it fast.
    </li>
    <li>
        The function <code>mandel()</code> is the simplest translation of
        the Mandelbrot fractal iteration into Python code that I could
        come up with.  The fact that Python has a built-in complex type
        makes it very easy to implement the standard Mandelbrot set
        algorithm.
    </li>
    <li>
        I have taken the maximum number of iterations
        inside <code>mandel()</code> to be 20, the same value used
        in the post I mentioned before.
        According to the very simple method used to time the application,
        it takes about 2 seconds on my computer to draw a
        simple picture.
        This is annoying slow.  Furthermore, by looking at the resulting
        picture, and trying out with different number of iterations
        in <code>mandel()</code>, it is clear that 20 iterations is not
        sufficient to adaquately represent the Mandelbrot set; this is
        especially noticeable when exploring smaller regions of the
        complex plane. A more realistic value
        might be to take 100 if not 1000 iterations which takes too long
        to be practical.
    </li>
    <li>
        Tkinter's canvas does not have a method to set the colour of
        individual pixels.  We can simulate such a method by drawing a line
        (for which there is a primitive method) of length 1.
    </li>
    <li>
        The screen vertical coordinates ("y") increase in values from the
        top towards the bottom, in opposite direction to the usual vertical
        coordinates in the complex plane.  While the picture produced is
        vertically symmetric about the x-axis, I nonetheless wrote the code
        so that the inversion of direction was properly handled.
    </li>
</ol>

<p>
    This basic application is not really useful as a tool for exploring
    the Mandelbrot set, as the region of the complex plane it displays
    is fixed.  However, it is useful to start with something simple
    like this as a first prototype.  Once we know it is working
    we can move on to a better second version.  So, let's write a fancier
    fractal viewer following the outline below:
</p>

<pre>
class Viewer(object):
    '''Base class viewer to display fractals'''

        # The viewer should be able to enlarge ("zoom in") various regions
        # of the complex plane.  I will implement this
        # using keyboard shortcuts.
        #
        self.parent.bind("+", self.zoom_in)
        self.parent.bind("-", self.zoom_out)
    def zoom_in(self, event):
    def zoom_out(self, event):
    def change_scale(self, scale):

        #  Since one might want to "zoom in" quickly in some regions,
        # and then be able to do finer scale adjustments,
        # I will use keyboard shortcuts to enable switching back
        # and forth between two possible zooming mode.
        # A better application might give the user more control
        # over the zooming scale.
        self.parent.bind("n", self.normal_zoom)
        self.parent.bind("b", self.bigger_zoom)
    def normal_zoom(self, event, scale=1):
    def bigger_zoom(self, event):

        # Set the maximum number of iterations via a keyboard-triggered event
        self.parent.bind("i", self.set_max_iter)
    def set_max_iter(self, event):

        # Like what is done with google maps and other
        # similar applications, we should be able to move the image
        # to look at various regions of interest in the complex plane.
        # I will implement this using mouse controls.
        self.parent.bind("<Button-1>", self.mouse_down)
        self.parent.bind("<Button1-Motion>", self.mouse_motion)
        self.parent.bind("<Button1-ButtonRelease>", self.mouse_up)
    def mouse_down(self, event):
    def mouse_motion(self, event):
    def mouse_up(self, event):

        # Presuming that "nice pictures" will be eventually produced,
        # and that it might be desired to reproduce them,
        # I will include some information about the region of the
        # complex plane currently displayed.
    def info(self):
        '''information about fractal location'''
</pre>



<ul>
    <li>
        Furthermore, while I plan to use proper profiling tools, I will nonetheless
        display some basic timing information as part of the GUI
        as a quick evaluation
        of the speed of the application.
    </li>
    <li>
        Finally, since I expect that both the function <code>mandel()</code>
        and the drawing method <code>draw_fractal</code> to be the
        speed-limiting factors, I will leave them out of the fractal viewer
        and work on them separately.  If it turns out that the profiling
        information obtained indicates otherwise, I will revisit
        this hypothesis.
    </li>
</ul>

<p>
    Here is a second prototype for our fractal viewer, having the features
    described above.
</p>

<pre>
''' viewer.py

Base class viewer for fractals.'''

import sys
if sys.version_info < (3,):
    import Tkinter as tk
    import tkSimpleDialog as tk_dialog
else:
    import tkinter as tk
    from tkinter import simpledialog as tk_dialog

class Viewer(object):
    '''Base class viewer to display fractals'''

    def __init__(self, parent, width=600, height=480,
                 min_x=-2.5, min_y=-1.5, max_x=1.):

        self.parent = parent
        self.canvas_width = width
        self.canvas_height = height

        # The following are drawing boundaries in the complex plane
        self.min_x = min_x
        self.min_y = min_y
        self.max_x = max_x
        self.calculate_pixel_size()
        self.max_y = self.min_y + self.canvas_height*self.pixel_size

        self.calculating = False
        self.nb_iterations = 20
        self.normal_zoom(None)

        self.canvas = tk.Canvas(parent, width=width, height=height)
        self.canvas.pack()
        self.status = tk.Label(self.parent, text="", bd=1, relief=tk.SUNKEN,
                               anchor=tk.W)
        self.status.pack(side=tk.BOTTOM, fill=tk.X)
        self.status2 = tk.Label(self.parent, text=self.info(), bd=1,
                                relief=tk.SUNKEN, anchor=tk.W)
        self.status2.pack(side=tk.BOTTOM, fill=tk.X)

        # We change the size of the image using the keyboard.
        self.parent.bind("+", self.zoom_in)
        self.parent.bind("-", self.zoom_out)
        self.parent.bind("n", self.normal_zoom)
        self.parent.bind("b", self.bigger_zoom)

        # Set the maximum number of iterations via a keyboard-triggered event
        self.parent.bind("i", self.set_max_iter)

        # We move the canvas using the mouse.
        self.translation_line = None
        self.parent.bind("<Button-1>", self.mouse_down)
        self.parent.bind("<Button1-Motion>", self.mouse_motion)
        self.parent.bind("<Button1-ButtonRelease>", self.mouse_up)

        self.draw_fractal()  # Needs to be implemented by subclass

    def info(self):
        '''information about fractal location'''
        return "Location: (%f, %f) to (%f, %f)" %(self.min_x, self.min_y,
                                                  self.max_x, self.max_y)

    def calculate_pixel_size(self):
        '''Calculates the size of a (square) pixel in complex plane
        coordinates based on the canvas_width.'''
        self.pixel_size = 1.*(self.max_x - self.min_x)/self.canvas_width
        return

    def mouse_down(self, event):
        '''records the x and y positions of the mouse when the left button
           is clicked.'''
        self.start_x = self.canvas.canvasx(event.x)
        self.start_y = self.canvas.canvasy(event.y)

    def mouse_motion(self, event):
        '''keep track of the mouse motion by drawing a line from its
           starting point to the current point.'''
        x = self.canvas.canvasx(event.x)
        y = self.canvas.canvasy(event.y)

        if (self.start_x != event.x)  and (self.start_y != event.y) :
            self.canvas.delete(self.translation_line)
            self.translation_line = self.canvas.create_line(
                                self.start_x, self.start_y, x, y, fill="orange")
            self.canvas.update_idletasks()

    def mouse_up(self, event):
        '''Moves the canvas based on the mouse motion'''
        dx = (self.start_x - event.x)*self.pixel_size
        dy = (self.start_y - event.y)*self.pixel_size
        self.min_x += dx
        self.max_x += dx
        # y-coordinate in complex plane run in opposite direction from
        # screen coordinates
        self.min_y -= dy
        self.max_y -= dy
        self.canvas.delete(self.translation_line)
        self.status.config(text="Moving the fractal.  Please wait.")
        self.status.update_idletasks()
        self.status2.config(text=self.info())
        self.status2.update_idletasks()
        self.draw_fractal()

    def normal_zoom(self, event, scale=1):
        '''Sets the zooming in/out scale to its normal value'''
        if scale==1:
            self.zoom_info = "[normal zoom]"
        else:
            self.zoom_info = "[faster zoom]"
        if event is not None:
            self.status.config(text=self.zoom_info)
            self.status.update_idletasks()
        self.zoom_in_scale = 0.1
        self.zoom_out_scale = -0.125

    def bigger_zoom(self, event):
        '''Increases the zooming in/out scale from its normal value'''
        self.normal_zoom(event, scale=3)
        self.zoom_in_scale = 0.3
        self.zoom_out_scale = -0.4

    def zoom_in(self, event):
        '''decreases the size of the region of the complex plane displayed'''
        if self.calculating:
            return
        self.status.config(text="Zooming in.  Please wait.")
        self.status.update_idletasks()
        self.change_scale(self.zoom_in_scale)

    def zoom_out(self, event):
        '''increases the size of the region of the complex plane displayed'''
        if self.calculating:
            return
        self.status.config(text="Zooming out.  Please wait.")
        self.status.update_idletasks()
        self.change_scale(self.zoom_out_scale)

    def change_scale(self, scale):
        '''changes the size of the region of the complex plane displayed and
           redraws'''
        if self.calculating:
            return
        dx = (self.max_x - self.min_x)*scale
        dy = (self.max_y - self.min_y)*scale
        self.min_x += dx
        self.max_x -= dx
        self.min_y += dy
        self.max_y -= dy
        self.calculate_pixel_size()
        self.draw_fractal()

    def set_max_iter(self, event):
        '''set maximum number of iterations'''
        i = tk_dialog.askinteger('title', 'prompt')
        if i is not None:
            self.nb_iterations = i
            self.status.config(text="Redrawing.  Please wait.")
            self.status.update_idletasks()
            self.draw_fractal()

    def draw_fractal(self):
        '''draws a fractal on the canvas'''
        raise NotImplementedError
</pre>
<p>
    We move the Mandelbrot set calculation in a separate file.
</p>

<pre>
# mandel1a.py

def mandel(c, max_iterations=20):
    '''determines if a point is in the Mandelbrot set based on deciding if,
       after a maximum allowed number of iterations, the absolute value of
       the resulting number is greater or equal to 2.'''
    z = 0
    for iter in range(0, max_iterations):
        z = z**2 + c
        if abs(z) >= 2:
            return False
    return abs(z) < 2
</pre>

<p>
    And, finally, we implement the missing functions for the viewer in
    a new main application.
</p>

<pre>
# viewer1.py

from mandel1a import mandel
from viewer import Viewer
import time

import sys
if sys.version_info < (3,):
    import Tkinter as tk
    range = xrange
else:
    import tkinter as tk

class FancyViewer(Viewer):
    '''Application to display fractals'''

    def draw_pixel(self, x, y):
        '''Simulates drawing a given pixel in black by drawing a black line
           of length equal to one pixel.'''
        self.canvas.create_line(x, y, x+1, y, fill="black")

    def draw_fractal(self):
        '''draws a fractal on the canvas'''
        self.calculating = True
        begin = time.time()
        # clear the canvas
        self.canvas.create_rectangle(0, 0, self.canvas_width,
                                    self.canvas_height, fill="white")
        for x in range(0, self.canvas_width):
            real = self.min_x + x*self.pixel_size
            for y in range(0, self.canvas_height):
                imag = self.min_y + y*self.pixel_size
                c = complex(real, imag)
                if mandel(c, self.nb_iterations):
                    self.draw_pixel(x, self.canvas_height - y)
        self.status.config(text="Time required = %.2f s  [%s iterations]  %s" %(
                                (time.time() - begin), self.nb_iterations,
                                                                self.zoom_info))
        self.status2.config(text=self.info())
        self.status2.update_idletasks()
        self.calculating = False

if __name__ == "__main__":
    root = tk.Tk()
    app = FancyViewer(root)
    root.mainloop()
</pre>

<p>
    Let me conclude with few black and white pictures obtained using
    this program, which, if you look at the time, highlight the need for
    something faster. First for 20 iterations, drawn in 4 seconds.
</p>

<img src="bw_20.png" />

<p>
    Then, for 100 interations - better image, but 7 seconds to draw...
</p>

<img src="bw_100.png" />

<p>
    Next post, we'll start profiling the application and make it faster.
</p>

</body>
</html>
